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ABSTRACT 


Highly maneuverable ocean going vehicles require quick response, control, and 
guidance to maintain accurate track keeping characteristics. Ocean vehicles, however, may 
experience significant lags in their inertial positional information that limit their reaction 
times. This thesis investigates the effects of these lags on guidance and control. The 
relationship of critical visibility versus the controller time constant and its effect on the 
stability of the guidance/control scheme is analyzed. Results are presented based on time 
domain and frequency response techniques using a dynamic model of the Naval 
Postgraduate School Autonomous Underwater Vehicle Il (NPS AUV II), for which a 


complete set of hydrodynamic coefficients and geometric properties is available. 
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I. INTRODUCTION 


A. BACKGROUND 

Marine vehicles that require high maneuverability, such as 
hydrofoils, require quick response, control and guidance to 
Memntain accurate track keeping characteristics. iiss 
accomplished through successful path merits navigation, 
guidance and autopilot design [Ref. 1]. 

Sufficient information is obtained from charted obstacles 
fmeaeechiec environment for smooth paths to be generated for the 
vehicle to follow [Ref. 2]. A certain level of feedback 1s 
provided through the use of sonar and acoustics in order to 
replan a path when uncharted obstacles are encountered or when 
the mission objectives are changed. Based on inertial 
positional information, the guidance law provides’~ the 
appropriate vehicle heading by suitable use of the vehicle 
effectors such as rudders, dive planes, and cross body 
thrusters. However, lags in obtaining and processing inertial 
MesttiOnal information can limit vehicle reaction time. In 
Pi@qaition, the guidance law must be as fast as possible in 
order to ensure accurate path keeping characteristics [Ref. 
3]. Therefore, the stability of the combined guidance/control 


scheme becomes an issue that needs to be analyzed. 


B. OBJECTIVE OF THIS THESIS 

This thesis investigates the effects of positional 
information time lags on guidance Vand) cone eu The 
relationship of the critical visibility (minimum vehicle 
lookahead distance) versus the controller time constant and 
its effect on the stability of the combined guidance/contres 
scheme 1s analyzed. Results are presented based on 
computations using time domain and frequency response 
techniques. All computations are performed for a Gyimemiaae 
model of the NPS AUV II [Refs. 4 and 5] for which a complete 
set of hydrodynamic coefficients and geometric properties 1s 


avallable. 


C. THESIS OUTLINE 

Chapter II presents the vehicle equations of motion in the 
horizontal plane, the equations that govern steering control, 
and the guidance scheme used in this research. 

Chapter III presents the stability analysis based on 
elgenvalue computation of the combined guidance and control 
scheme developed in Chapter II. Results are presented based 
on the stability curves developed for first, second and third 
order approximations for time lag in the commanded vehicle 
heading in the control law. 

Chapter IV presents an exact computation of the stability 
curves based on frequency response methods, which utilizes the 


Nyquist criterion for stabaiae7 


II. EQUATIONS OF MOTION, CONTROLLER DESIGN AND GUIDANCE SCHEME 


A. INTRODUCTION 

The equations of motion that describe maneuvering of a 
marine vehicle in the horizontal plane are presented in this 
chapter. A linear full state steering feedback control law is 
designed based on the three equations in sway, yaw, and rate 
of change of heading angle. The guidance scheme is based on 
pure pursuit guidance for path following along straight line 


segments. 


B. EQUATIONS OF MOTION 

fe -eGleeimgmour AQttemnlon tO Motion in the horizontal 
plane (steering control), the mathematical model consists of 
the nonlinear sway and yaw equations of motion only. Witha 
moving coordinate frame fixed at the vehicle’s geometric 


center, Newton’s equations of motion are 


MGQVer tat Xan de—/ ) = Y [Ze 


Tacit Valtie env r —N (2.2) 


where v and r are the relative sway and yaw velocities of the 
moving vehicle with respect to the water; mis the mass of the 
vehicle; x, y, are the respective lateral and longitudinal 
positions of the center of gravity; and Y,N represent the 
total excitation sway force and yaw moment, respectively. 
These forces can be expressed as the sum of quadratic drag 
terms and first order memoryless polynomials inv andr, which 
represent the added mass and damping due to the vehicle’s 
motion through the water. In this way, Y and N can be 


represented by 


Y = it. (¥,v+¥_ur) +Eey uv 


- 2 ae 5) LED pags Beyure 


tail me | Vireo 


N= _ ORE Aci vue) +LONuv 


_ Bc h g) tbr) edt+ Penu7s 


taiz PY Jv+Er| 


where @ is the vehicle length, and Y,, N, represent partial 
derivatives of Y and N with respect to a, ©), 25.0 =suem 
coefficient, and 6 1S the ruddensanqser 

Equations (2.1) and (2.2) can be nondimensionalized wie 
respect to the constant forward speed u, and the vehicle 


length @; with the dimensionless time variable being tu/@. 


The cross flow integral drag terms become important for 
hovering operations or low speed manuevering, whereas at 
relatively high speeds and low angles of attack (with respect 
to the water), the steering response 1s predominantly linear. 

The model becomes complete with the addition of the 


expressions for the vehicle yaw rate 


Ww=r (2535) 


and the inertial position rate (horizontal plane) 


y=usiny+vcoswp (2a) 


where Wis the heading angle as shown in Figure 1. 
Pemabteons (2.1), °{2.2), and (2.3) can be written as a set 


of three coupled linear differential equations of the form 


w=r (AR Ses, 
v=a,,uvt+a,,ur+b,u*d (2 .9b) 
aa tiv. tt 11-0 esc) 


Figure 1. Vehicle Geometry and Definition of Symbols. 


where the coefficients a,; and b; are functions of the vehicle 


hydrodynamic and geometric properties, wis the heading angle 


Pima O 1S the rudder angle. 


Mice linearized form of Equation (2.4) 


yeup+v 


Hence, the final set of linearized equations of motion 


for steering control are 


Ww=r 


; 2 
V=a,,uvta,,ur+b,u*d 


_— 2 


yeupty 


at any nominal u. 
This system of equations becomes 


Sencroller design. 


ie Ieyelstaye 


apres nominal ai is 


used 


ea) 


7D) 


wie) 


viel), 


the 


Ce CONTROLLER DESIGN 
Equations (2.7a), (2.7b), and (2.7c) govern the steering 
control of the model used in this teseqmaa 


These equations can be represented in the form 


X=Ax+b (2.8) 


where the state vector equation 1s 


Pe [Weve (2.9) 


Linear full state feedback is introduced in Ghee 


S=k, (W-.) tk, v+k,r (2 Ge 


where yw. 1s the commanded heading angle, and the gains k,, kj, 
k, are computed by pole placement such that the chosed@ileag.s 
system of equations (2.8), (2.10) has the desired dynamics. 
The existence of the difference (w-yw.) in the control law 
(2.9) has the effect of trying to point the vehicle’s 
longitudinal axis towards the desired heading. 

If the desired characteristic equation is selected so that 


its time constant is t. dimensionless seconds, its general 


form becomes 


(woah 1 Wea ey Geel Rei Cae) 


where 





3 ee ik 
ere 1 on 5 2 rs 
c Ge eA 
The controller gains are computed from 
a 
—————— (2.12a) 


(b, a,,—b,a,,) uk, + (b,a,,—b,a,,) uk, =@,+b,u?k, = (4. 12b) 
b,u*k,+b,u°k,=-@,- (a); +422) u (2.12c) 
Due to the explicit dependence on u, the gains k,, k,, k, are 


continuously updated for different nominal forward speeds. 


1: GUIDANCE SCHEME 
The guidance scheme used in this research 1S an 


Orientation based command scheme. CMe eomecehneme. che 


commanded vehicle heading angle w. 1S a function of the line 
of sight angle 6 between the actual vehicle position and a 
reference position on the nominal path which is located at a 
constant lookahead or visibility distance d ahead Gi 
vennoe as shown in Figure 1. 


This line of sight angle 1S detimedme 


tang=-— (27 


AIS 


The autopilot is then called upon to deliver the commanded 

heading angle w.. The simplest orientation guidance law is 
pure pursuit guidance where the commanded heading angle yw, in 
the control law (2.9) equals the line of sight angle oO in 


Cobos 


Then yw. 1s defined by 
= -arctan (2. 


For relatively small angles 


y= -arctan4 = — (200s 


10 


This results in a three state heading autopilot design 
which can exhibit very robust characteristics. However, the 
stability of this scheme may become an issue when transiting 
along straight line segments since the commanded heading angle 
is a function of the vehicle response, l1.e., autopilot 
response is limited by time lags in vehicle dynamics. We will 


next examine the stability of this scheme. 


iGac 


III. EIGENVALUE ANALYSIS 


A. INTRODUCTION 

In this chapter, we present the stability analysis, based 
on eigenvalue computation, of the combined guidance and 
control scheme developed in Chapter II. A brief presentation 
of the stability curve is developed first for the case of zero 
time lag. Incorporation of a nonzero time lag eae 
positional information formally results in an infinite state 
space system. This is truncated into first, second, andi@iagms 
order approximations. The resulting eigenvalue problems are 


then analyzed with emphasis on stability curves. 
B. STABILITY CONSIDERATIONS AND TIME LAG 
1. Stability Considerations 
In pure pursuit guidance, where the commanded heading 
angle w. in the control law equals the line of sight angle Oo, 


the trivial equilibrium state which corresponds to straight 


line motion is characterized by 


W=v=r=y=0 (3 


Linearization of the state equations in the vicinity of (3.1) 


12 


produces the linear system 
X=AX, (332 } 
where the complete state vector 1s 


a= i ie, Vile (S537) 


Local stability properties of (3.1) are then established by 
the eigenvalues of A. The characteristic equation of (3.2) is 


meeeuartic of the form 


44+BA2+CA2+DA+E=0, (3.4) 


where the coefficients B, C, D, E are functions of the vehicle 
properties, the lookahead distance d, and the controller gains 
i[ess, K,- A pair of complex conjugate roots of (3.4) crosses 
the imaginary axis when the following conditions are met 


Meer. 6]. 


BCD-B*E-D=0, (3.5a) 


Pyo (3.5b) 
B 


13 


These two conditions (3.5a) and (3.565)) Crane ee 


a,d*+a,d+a,=0, (3 .6a) 


d> Db, a,,-b,4,, -D, (3 6e 
b,a,,—D,a,, 
where 
a,=@,a,-a, (3 .7a) 
_ £450,725) (Drage “Drerz by) __ yay _ gay 
° D, aya (b,a,,-b.a,,)u 1 (ae 


a,=— (D, 4,,—b,4,,~b,) [b, 4, + (b,a,,—b,a,,~b,) ul a, (3.7c) 
(b,a,,~b,a,,) *u 


The conditions (3.5a) and (3.5b) determine the critical 
visibility ds. for stability. For d > d.i., the equa iam 
state (3.1) is stable which means that the control law will 
drive and keep the vehicle onto the straight line path. For 


Sy cd HO! the equilibrium state loses its stability and the 


Sr pee 


vehicle response becomes oscillatory as a result of the pair 


14 


of complex conjugate eigenvalues with positive real parts. 

fro (lies oLmeemC © Cittleal visibility versus the 
controller time constant t. for zero time lag are presented in 
Figure 2. As expected, the higher values of t, (1.e.,softer 
controller) require higher lookahead distances d for path 
accuracy. This agrees with existing guidelines that the 
guidance law must be sufficiently slower than the controller 
in order for the dynamics of the two not to interfere with 
each other and, therefore, guarantee stability of the combined 
scheme. Very high values of d correspond to a very slow 
guidance law with a loss in speed of response and path 


pecuracy [Ref. 1]. 


2. Time Lag 
In ocean vehicles, all states needed in the control 
law are readily available at the required rate, with the 
possible exception of the positional information y (Figure 1). 
The latter may require a significant data analysis and 
reduction of sonar returns and inertial navigational 
information. As a result, it is likely that a time lag will 
Peat in the positional Porat Lon y which is used in the 
guidance law. 
With the introduction of a time lag T (sec), the commanded 


rudder angle yw. in the pursuit guidance control law becomes 


LS) 


2.0 


0 0.5 


Figure 2. 
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Distance d versus t, for 


Vv. = -arctanVie (35-3°) 


and for relatively small angles 


fi eM (3.9) 
mee resultant rudder control law (2.9) becomes 
b=k, wrk, LED ek vekyr eeun 


After some algebra, the resultant linearized equations of 


mememon (2.7a), (2.7b), (2.7c), (2.7d) become 


W=r (oo, 1 la) 
; 2 2 2 2 ky 
vou Kt (a,,u+b,u°k,) v+(a,,u+b,u?k,) r+b,u ay (e-T) Va od lisy, 


k 
f=b,u*k,p+(a,,u+b,u7k,) v+ (a,,u+b,u7k,) r+b,u?——y( a oC) 


7 


yr=uytv (3.1 ee 


The Taylor expansion for the ternage wa 


2 3 
vile ie Y-TY+ —Y-Ie - (3. dee 


We can now examine the stability of the previously developed 
Crore 1010) 1 scheme FOr Fil eete second and thine order 


approximations for time lag in the contra ler = 


C. FIRST ORDER APPROXIMATION FOR TIME LAG 


For a first order approximation for time lag T 
VV 6-f) ya (3 aan 


where 


y=uyt+v (30s) 


After some algebra, the resultant linearized equations of 


motion (3.lla), (3.11b), (3 Lic} (Gye weeeen 


18 


er 


(ay ibscu) 
‘ 2 cage b 2k -b 2X1 + ( +b 2k.) +b 2 Ky 
v= (b,u Ie ee | qT) + (aut pK), “d ONS NS eae ie Pa eel ob) 
, 2 3 ky 2 2 k, 2 b.u2 ky 
f= (b,u?k,~b,u* — T) p+ (a,,u+b,u*k, -b,u*— T) v+ (a,,ut+b,u°k,) r+b,u"-F¥(3 15) 
y=uytv (3 i Sah) 
The state space form of the equations of motion is 
xX=AX (25.16) 
mamm@ieter ix corm, (3.15a), (3.15b), (3.15c), (3.15d) become 
y 0 0 ft 0 ¥ 
‘ 2 3K 2 2 Ky 2 2 Ky 
v Contre Du qi (a,,u+b,u?k,-b,u qr) (a7,U+D,u-k,) (blu Fil V (3.17) 
18 (bju?k, ~b,u? 3 7 (a,,u+b,u?k,-b,u?=27) (a,,u+b,u*k;) (bju?<3) 36 
ly} : " . . Ly 


The standard eigenvalue problem (3.17) is solved to 
determine numerically the critical d versus t,. curve for a 


Given value of T. Appendix A presents the program used for 


SEY, 


the first order approximation f£Om mame eae 

Figure 3 shows the resulting stability curves for the 
first order case. It shows the relation between the critical 
visibility d and the controller time constant t. for Vali@eewem 
time lag of 0.0, 0.5, 1.0, 1.5 amad™ 27 0Ne cect 

As expected, the stability curve shifts to the right with 
increasing values of time lag. For the same time constant 7; 
the critical visibility distance increases with corresponding 


increases in the amount of time lag. 


D. SECOND ORDER APPROXIMATION FOR TIME LAG 


For a second order approximation for time lag T, we have 


2 
y(e-T)* y-Ty+2y See 
where 
y=uyrtv, yr=urtv (3215 t (3 as) 


Following some algebra, the resultant linearized equations of 


motion (3.lla), (3.11b), (3.11¢c)) (30 ieee eee 


Zz) 


te 


Figure 3. 





First Order Approximation For Time Lag: 


Critical Visibility Distance d versus t, for 


Time Lags T = 0.0, 0.5, 1.0, 1.5, 2.0 sec. 


BAe 


=r (3 .20a) 


V=A,,W+A,,V+A,,0°+AsV (3.205) 
D=A,, WtAl.V Ao eee (3.2068 
y=uytv (3 .20d) 
where 
k 
b,u?k, -b,u? = T 
ay = iE ' 
1-b,u?—=T? 
1" 2d 


k 
aan oaks “bu? iG 


k 
1-bu?* >? 


kK, 
a, urbyudka baa oF hs 


k 
1-b,u* 5-7? 


Ze 


(b, u?k, -b,u? af T) 
k, Ke SS Se 
if : d 2d oe 
(1 bu Bae ) 


k 
2k, -b,u? = T 
a. a ie, | oe (2a am erie Ks late F ) 
a pa eros | k, lop at + (b,u* —T po , 


2d = BSS 
ah T?2 


k 
(a,,u+b,u*k, tb,ue 5-7?) 


k 
= 2 al 2 
(Gl eke sar ) 


k 
A;; = (a,,u+b,u*k,+b,u*——T*) + (bau? 5-7?) 


As with the first order case, the state space equations of 


Mmeebron reduce to 


Tei 


X=AX (3 . ie 


In matrix form, (3.20a), (3.206), (3.20c), (3.200) seee7— 
y 0 0 a 0 Y 
V Az, Az2 Ap; As, VY 
= (3 2a 
r Az, A32 A333 Asy r 
MK u att 0 0 , 





Simular to the first order case, the standard eigenvalue 
problem (3.21) 1s solved using Fortran programming and MATLAB 
in order to develop the stability characteristics of the 
second order case. Appendix B presents the program used to 
develop the stability characteristics for this case. 

Figure 4 shows the resulting stability curves for this 
case. The results of the critical visibility d versus the 
controller time constant t,. for time lags of 0.0, 0 3Syye 
1.5 and 2.0 secs are shown. 

As with the first order case, the higher values (Glue 
require higher lookahead distances d for path accuracy. High 
values of d correspond to a slow guidance law with a loss in 
speed of response and path accuracy, and the stability curve 


shifts to the right with increasing values of time lag. 
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Figure 4. 





Second Order Approximation For Time Lag: 


Critical Visibility Distance d versus t, for 


hime wbags 1 = 0-20, 0.5, 1.0, 1.5, 2.0 sec. 


A 


However, due to being more accurate than the first order case, 
the second order approximation provides tighter stabil 


curves for the same time lags. 


E. THIRD ORDER APPROXIMATION FOR TIME LAG 


For a third order approximation £6r timeoelaaee 


2 ) 
Usa) = y-Ty+ = y-2V (3.208 
where 
yruy+v, yrturtyv, yrurt+yv (3.14), (3.19), Ga 


After a considerable amount of algebra and following a simular 
procedure as with both the previous two cases, the resultant 


linearized equations of motion reduce to the state space form 


BX=AX (3 245 


Cs 


X= Bees 
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The resultant linearized equations of motion 


meio), (3.lle), (3.l1ld) become 


B,,2+B,,V, = A, W+A,,V,+A,,V,+A,,°+A, ov 


Be30+Be5Vz = As, WtAc, V1 +As53 V2 +Acgr tAcc¥ 


y=uyrv, 
Vv, = Vv; 
where 

k 

B,, = bu’ T° 
k 

B,. = bu? T 

2 3 Ky 
Een ao) ibs ai 


Prat 


(e. 


TP a)e, 


25a) 


Eye) 


2251C) 


oe) 


.25e) 


k 
A,, = @),U+b)U 2 een oe 


k 
Az4 = b,u?— 


k 
alee = bu? 77-1 


and 


k 
Ac, = b,u*-b,u*—-T 


k 
Ac, = a,,U bynes ~b,u*— iy 
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k 
A,, = a,,u+b,u*k,+b,u?—=T’ 


2a 
k, 
As, = bu*— 
k, 
A,, = bu*>.T* 
Mem i< rom, eee na)> (5.255), (3.25c), (3.25d), (3.25e) 
become 
1 Se Re y 0 Ti ae y 
oS Cmte vy, 0 oOo oOo oO 1 Vi 
Geo By. 150 B,. pW GAs Ase? GA. “Ae “Aj. r}{ (3.26) 
0 0 0 al 0 y u i 0 0 0 y 
0 O 38, 0 Bss V2 As, As, As; Ase Ass Ve 


Simular to both the first order and second order cases, 
the generalized eigenvalue problem (3.21) is solved using 
Fortran programming and MATLAB in order to develop the 
stability characteristics for the third order case. Appendix 
C presents the program used to develop the stability 
@earacteristics for this case. 


Figure 5 shows the resulting stability curves for this 
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Figure 5. 
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Third Order Approximation For Time Lag: 
Critical Visibility Distance d versus t, for 


Time Lags T = 0.0, 0.5, 1.0, 1.5, 220eece 
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Peewee jes rectitrs Of the Critical visibility d versus the 
Bontroller time constant t. for time lags of 0.0, 0.5, 1.0, 
1.5 and 2.0 seconds are shown. As with the first and second 
order cases, the higher values of t. require higher lookahead 
Mestances dg for path accuracy, high values of d correspond to 
a slow guidance law with a loss in speed of response and path 
accuracy, and the stability curve shifts to the right with 
increasing values of time lag. However,: due to being more 
accurate than either of the other two cases, the third order 
approximation provides the tightest stability curves for the 


same time lags. 


F. COMPARISON OF RESULTS 

Results of the critical visibility d versus the controller 
mime) Constant t. for first, second and third order 
approximations for time lag have been obtained. Figures 3, 4 
and 5 presented the stability curves for first, second and 
merc Order approximations for time lags of 0.0, 0.5, 1.0, 1.5 
and 2.0 secs. These figures have shown that stability of the 
control scheme decreases with increasing time lag, and for the 
Same time constant, the critical visibility distance increases 
with increasing time lag. 

Figures 6, 7, 8 and 9 present a comparison of first, 
second and third order approximations for each time lag 
considered. These figures show that for each time lag, the 


third order model presents the largest region of stability for 
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Figure 6. 


_istord |  , 7 
2nd Ord 
3rd Ord 


“ Prete ere etre eee ery eee er vere te eee 


Critical Visibility Distance d versus t, for 
First, Second and Third Order Approximation 


Cases and Time Lag T = 0.5 sec. 
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Figure 7. Critical Visibility Distance d versus t,. for 
First, Second and Third Order Approximation 


Cases and Time Lag T = 1.0 sec. 
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Figure 8. 
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Figure 9. Critical Visibility Distance d versus t, for 
First, Second and Third Order Approximation 


Cases and Time Lag T = 2.0 sec. 


a. 


vehicle straight line path accuracy. In addition, £o0r ge 
same time constant, the critical visibility is least fOr wie 
third order case. This demonstrates that the minimum 
lookahead distance required for straight line motion stability 
decreases with controller time lag accuracy. 

It can be seen that the first order approximation Geum 
most conservative for design purposes since it predicts the 
highest required minimum visibility distance. As expected, 
the differences among the three approximations are more 
pronounced as the time lag increases and as the controller 
time constant becomes smaller; 1.4€., tighter Gonuueme 

Figures 10, 11, 12 and 13 present the differences in 
critical visibilities among first, second and third Wome 
approximations for each individual time lag considered. A 
comparison of these figures demonstrates that the difference 
in critical visibility between respective orders (time lag 
models) decreases with decreasing time lag. This indicates 
that for low time lags, controller time lag accuracy is not as 
crucial for maintaining vehicle straight line path accuracy as 
that of higher time lags. The greater the time lag, the 
greater the difference in critical visibilities among time lag 
models, and the more crucial is the controller time lag 
accuracy for stabrllry= 

In general, if the guidance law is slow enough to allow 
sufficient time for the controller to react, stability of 


Straight line motion is guaranteed [Ref. 1]. 
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Figure 10. Difference of Critical Visibility Distances at 
Time Lags of T = 0.5 and 0.0 sec versus t, for 


First, Second and Third Order Cases. 
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Figure 11. Difference of Critical Visibility Distances at 
Time Lags of T = 1.0 and 0.0 sec versus t, for 


First, Second and Third Order Cases. 
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Figure 12. Difference of Critical Visibility Distances at 
Time Lags of T = 1.5 and 0.0 sec versus t, for 


First, Second and Third Order Cases. 
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Figure 13. Difference of Critical Visibility Distances at 
Time Lags of T = 2.0 and 0.0 sec versus t, for 


First, Second and Third Order Cases. 
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G. NORMALIZATION OF RESULTS 

In an attempt to characterize the stability analysis 
Mmesults for all values of time lag, we proceed as follows: If 
Mme vehicle side slip velocity v is very small so that it can 
memmecaglected, the equations of motion (2.7a), (2.7b), (2.7c), 


mma (2./d), take the form 


a 
li 
Hi 


-— 2 
ee eM 


y=uy 


m—memee vy = 0. The control law is 


d=k, (W-W.) +k,r, 


where the first order approximation for the commanded heading 


angle is 


Based on the above equations, we can get the characteristic 


Al 


equation of the system as 


s?-(a,,u+b,k,u?) s?-(b,k,u?-b,k,u3?T+) s-b,k,u3 


A = 00 


= 
d 


Applying Routh’s criterion to this cubic equation, Wwe siaaae 


that loss of stability occurs when 
(a,,u+b,k,u2) (b,k,u?-b,k,u3T+) = -b,k,u324, 


d d 


from which we can get the critical visibility distance 


a ee 5) (3.2 


Equation (3.27) can be written as 





y. (3.288 


where d, 1S the critical visibility distance for T = Uj@=amaueee 
the critical visibility distance for a nenZzemeusas Ligh 


dimensionless quantities, equation (3.28) is 
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= Goa 9) 





Mere 1 = 2.0 £L/sec 1S the vehicle speed, and @ = 7.3 ft is 
the vehicle length. 

A comparison among the approximate expression (3.29) and 
Me@ewtirst, second, and third order approximations is shown in 
figures 14, 15 and 16. The agreement is better for the first 
order approximation, as expected, and also for the higher 
controller time constant t.. This is because a high t, results 


in soft vehicle response with negligible amounts of side slip. 
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IV. FREQUENCY RESPONSE ANALYSIS 


A. INTRODUCTION 

The previous analysis using Taylor expansions of y(t-T) 
breaks down for approximations beyond third order. The reason 
Bor this is that for higher than third order expansions, the 
matrix B in the generalized eigenvalue problem (3.23) becomes 
Singular. Therefore, 1t was felt that a different technique 
should be sought in order to obtain an exact computation of 
the stability curves and also to check the validity of our 
calculations. In this chapter, we present a technique based 
on frequency eens methods, which utilizes the Nyquist 
M@eeeerion for stability. The resulting equation is solved 


uSing Newton’s iteration. 


B. NYQUIST CRITERION 
From our previously developed model, the linearized 


horizontal equations of motion are 


wH=r (2.7a) 


v=a,,uvta,,ur+b,u*d (27715) 


4'7 


r=a,,uvt+a,,ur+b,u7s (2 . 7am 
y=uyt+v (2.7a8 

at any nominal u. The control law is 
5=k, (W-W,) +k, vt+k,xr (2 oF 


where the commanded heading angle is 


ge SRD 346 
Wo ae (3 399 


The Laplace transform of equation (3.9) is 


eis (4.1) 


and the resultant rudder control law (2.10) becomes 


k 
b=kp+k,v+k,r+—-ye™ (4.2) 


Following some algebra, the resultant linearized equations of 
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fmeteme@n (2.7a),e(2275), (2e7c), (2./d) become 


(43a) 


<= 
it 
NX 


k kK 
v= (b,u?k, -byu? 27) yr (a,,u+b, u?k, ~b,u?— T) v+ (a,,u+b,u7k,) r+b,u’—ye™ (4.3b) 


b= (byu?k,-byu? 27) yr (4,,u+bju?k,-b,u? 27) v+ (a,,u+b,u7k;) r+b,u* = ye-™ (4.36) 


y=uytv ( 73) 


fmetetons (4.34a)7" (4.35), (4.3c), (493d) reduce to the state 


Space form 


= Qe) 
The matrix form becomes: 

¥ 0 0 1 0 ¥ 

V Corde la ued, uk) (a,,U+D,u-K,) (b,u2 Xt. 9-t5) 

d 
a (4.4) 
a (ue eas utpoU- i) “(aoju+b,u*k,) (byu? he -T8 es 
> | . 0 0 pu | 


AQ 


The characteristic equation of the form [{A-Is] = 0 began 


after a considerable amount of algebra 


S‘+A,s'+A,s7+A,S+(B)s@ 7 ese, es an (4.5) 


where 


A, = ~(4,,+42,) u -(b,k,+b,k,) u* 


A, = (4, 1452—4,24>,) uae (a,,5,-a,,5,) u*k,+ (a,,b, =aipoe Bes, -b,us 


A, = -(a,,b,-a,,D,) u°K, De *"=(a>, bya. mes 
B= =p 
B= bank 


and 


Ale 
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The characteristic equation of the system is written as 


1+kK G(s) =0 Ee 


where 


Dates Bisse ie 
S(s?+A,s*+A,st+A,) 


Memrne Eranster function, and we have denoted 


G9) (4.8) 


With s=jW, the phase angle 0 1S represented by 


b = £G(jo) (29) 


meet LOollows that 


ob = pie. ®) — LG@) + AB Oo. 7 O48.) - VO. w- AL 7+ A.) 
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For stability, the Nyquist criterion states (Figure 17) 


at the solution of the equarzven 


d(w) = -7% (4.11) 


which is phase cross over frequency @ = @,, the gain margin 


must be equal to l 


|KG(jw,)| =1 (4 ae 


where K =D = 1/d. It follows that equation (4.11) becomes 


3 
@ -@)) +A, W 
-w,T-= +tan™? (—+—+_) -tan"?(———4—=) =-n (4.13) 
By BW} -A,@,+A, 
Rearranging, we get 
Ww W>+A,W 
tan™* (-—-—>.) -tan™* ( = —— -T+—+0,T 
By B, Wi -A, W +A, 


Taking the tangent of both sides and using the identity 


tan (x) -tan(y) 


tan (x-Y) Sean ees 


BV 


v 





GH(jw) | 


Figure 17. Example of the Nyquist Stability Criterion For 


Three Different Values of Gain [Ref. 7]. 


op 


and rearranging following some algebra, the result becomes 


B,w, (A,-A,@}) -(By-B,@{) (-w]+A,0,) | : 
(B,-B,w{) (-A,{+A,) +B,o, (-@{+A,0,) tan(w,T) (4.14) 


From equation (4.12), 1 tolleyvssenee 


2 ' 

eo |-B,w1+B,jwo,+B,| =] 
aoe 

e ©, |-J01-A,0}+A,j0, +A, | 


Solving for ad, we Gee 


7) ae: 
Pe (55-2501) °teso (4.15) 


Ww, (A,-A,w{) 2+ (A,@,-@})? 


With a time lag T = 0, we get from equation (4.14) 


(B,-B,W{) (A,w{+A,) + B,w,(-W}+A,0,) = 0 


Rearranging, we get 


(B,A,-B,) w} ° (B,A,-B,A, -B,A,) 0; + BA, = 0 (4. Ge 
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Equation (4.16) can be solved exactly for @,, and then d 
ea be computed from (4.15). This zero time lag calculation 
was performed in order to validate the results of d versus t, 
/meeeche zero time lag solution of equation (3.6). The results 


using the two different techniques were found to be identical. 


C. ALGORITHM DESCRIPTION 


Introducing the terms 
@, = BA,-B, 
a, = B,A,-B,A,-B,A, 


=A, 


and 
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Then equation (4.12) becomes 


B,0°+B,0°+B,0 a ils 


@,W°+a4,07+4, tan (wT) 


Rearranging, we get 


(B,@°+f8,0°+B,0) tan (@7) = Go@- +50 "eee 


Fquation (4.18) is now tavine Stour 


Following Newton’s iteration for W@W, we have 


0,2 0,, - a 
f° (W,.4) 
From equatwven, (4.7008 
f(w) = (P,o0°+B,0°+B,o) sin (wT) 


+ (a,0°+a,07+a,) COS (WT) 
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(4.17) 


(4.18) 


(4.19) 


(4.20) 


(4.208) 


and 


=e (wo) 


(5B, 04+3B,w*+B,) sin (wT) 


+ 


(B.w°+B,w2+B,0) Tcos (wT) 
pee ae 2s 
+ (4@,w0°+2a,w) cos (wT) 


- (a4@,w*t+a,w7+a,) Tsin (wT) 
Ale 2 3 


Fortran programming and MATLAB were used in order to 
develop the stability characteristics for the Newton’s 
iteration method. Appendix D presents the program used to 
Gevelop the stability characteristics for this case. 

Figure 18 shows the resulting stability curves for the 
Newton’s iteration method. ijewesscmltis of Ehe critical 
feeb ility ad versus the controller time constant t. for time 
meemperor 0.0, 0.5, 1.0, 1.5, and 2.0 aes are shown. 

As with the first, second and third order approximation 
cases, the higher values of t. require higher lookahead 
distances d for path accuracy, high values of d correspond to 
a slow guidance law with loss in speed of response and path 
accuracy, and the stability curves shift to the right with 
increasing values of time lag. However, since the Newton’s 
iteration method presents the exact solution of the frequency 


response to equation (4.19), the resulting stability curves 
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Figure 18. 
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Newton’s Iteration Method: Critical Visibility 
Distance d versus t, for Time Lags T = 0.0, 


0.5, 1.0, 2.0 secs. 
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represent the precise stability characteristics for these time 


lags. 


D. RESULTS AND DISCUSSION 

Pe he csOmeriererMElcal Vvisibl lity d versus the controller 
time constant t. for the Newton’s iteration method has been 
demonstrated. Figure 18 presented the precise stability 
mitreves for time lags of 0.0, 0.5, 1.0, 1.5 and 2.0 secs. 

This figure, as well as those for the first, second and 
third order approximation cases, demonstrates that the 
stability of the control scheme decreases with increase in 
time lag and for the same time constant, critical visibility 
increases with increase in time lag. 

Figures 19, 20, 21, 22 and 23 present a comparison of 
generated stability curves among the Newton’s iteration method 
and those for the first, second and third order approximation 
Mmieo tor time lags of 0.0, 0.5, 1.0, 1.5 and 2.0 secs. It 
can be seen that there is barely a slight difference in the 
stability curves for the Newton’s iteration method case and 
the third order approximation case. This indicates that a 
mora Order approximation for time lag presents the best model 
MemmeccCOUNtLing For a positional information time lag in the 


melicle control law. 
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Figure 19. Critical Visibility Distance d versus t, for 
Newton’s Iteration Method and First, Second and 


Third Order Approximations; T = 0.0 sec. 
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Figure 20. 
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Critical Visibility Distance d versus t, for 
Newton’s Iteration Method and First, Second and 


Third Order Approximations; T = 0.5 sec. 
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Figure 21. Critical Visibility Distance d versus t, for 
Newton’s Iteration Method and First, Second and 


Third Order Approximations; T = 1.0 sec. 
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Figure 22. Critical Visibility Distance d versus t, for 
Newton’s Iteration Method and First, Second and 


Third Order Approximations; T = 1.5 sec. 
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Figure 23. Critical Visibility Distance d versus t, for 
Newton’s Iteration Method and First, Second and 


Third Order Approximations; T = 2.0 sec. 
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CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

A methodology for considering positional information time 
lags in the control law for vehicle maneuvering in the 
horizontal plane has been presented. The relationship of the 
Mame rcal Vasibility versus the controller time constant for 
both zero and non-zero time lags has been established. 

Time lags have been approximated by first, second and 
third order Taylor expansions in the commanded vehicle heading 
Mmmeme Gomntrol law. A comparison of the resulting stability 
SUrves has demonstrated that the third order approximation 
presents the greatest stability and the least critical 
visibility for a given time lag. It was found that the first 
order approximation was the most conservative for controller 
design purposes since it predicted the highest required 
Sem@etcal visibility distance. In addition, the differences 
among the three approximations became more pronounced with 
increasing time lag and decreasing controller time constant 
memohter control). 

Further analysis was conducted to characterize the 
stability results for any value of time lag using a 
normalization procedure and the first order approximation for 


time lag in the commanded vehicle heading in the control law. 
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A comparison of the normalized stability curves for Gach 
lag showed that the greatest differences occur with increasing 
order and decreasing controller time constant. This is a@W@euee 
a tight (quick) vehicle response with Significant side slip at 
lower time constants. 

Frequency response techniques based on the Nyquist 
Stability criterion were used to produce the exact stability 
curves for comparison with the previously generated Taylor 
expansion stability curves. These results virtually matched 
those obtained for the third order approximation time lag case 
and validated our calculations and stability curves generated 
for the case of zero time lag. These results indicateameaem 
the third order approximation best represents a positional 
information time lag in the control law, however for design 
purposes, the first order approximation provides the most 


conservative approach. 


B. RECOMMENDATIONS 
Some recommendations for further research are as follows: 

‘ Experimental verification using the NPS AUV II. 

a Consideration of time lag in the simulation of an Inertial 
Navigation System required for positional updates. 

= Analysis of positional information time lags in the motion 
Stability in the vertical plane, along curved paths, or 


with respect to other guidance schemes. 
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APPENDIX A 


PROGRAM THESIS1.FOR (Time Delay~-lst Order Approx) 
BIFURCATION ANALYSIS 


IMPLICIT DOUBLE PRECISION (A-H,0O-Z) 

DOUBLE PRECISION K1,K2,K3,L,NR,NV,NDRS,NDRB,IZ,MASS, 
& NRDOT,NVDOT 

DIMENSION A(4,4),FV1(4),IV1(4),222(4,4),WR(4) ,WI(4) 


OPEN (11,FILE='BIF1.RES’ ,STATUS=‘NEW’ ) 
OPEN (12,FILE='BIF2.RES’ ,STATUS=‘ NEW’ ) 
OPEN (13,FILE=’BIF3.RES’ ,STATUS=‘ NEW’ ) 


Vehicle Parameters: 
WEIGHT=435.0 


IZ =45.0 

L =7.3 

RHO =1.94 

G =32.2 

XG =0.0104 
MASS =WEIGHT/G 


YRDOT =-0.00178*0.5*RHO*L**4 
YVDOT =-0.03430*0.5*RHO*L* * 3 
YR =+0.01187*0.5*RHO*L** 3 
YV =-0.03896*0.5*RHO*L**2 
YDRS =+0.02345*0.5*RHO*L**2 
YDRB =+0.02345*0.5*RHO*L**2 
NRDOT =-0.00047*0.5*RHO*L**5 
NVDOT =-0.00178*0.5*RHO*L**4 


NR =-0.01022*0.5*RHO*L**4 
NV =-0.00769*0.5*RHO*L* * 3 
NDRS =-0.337*0.02345*0.5*RHO*L** 3 


NDRB #+0.283*0.02345*0.5*RHO*L**3 


WRITE (*,1001) 
READ (*,*) 

WRITE (*,1002) 
READ (*,*) 

XDMIN=XDMIN*L 
XDMAX=XDMAX *L 
WRITE (*,1003) 
READ (*,*) U 


TCMIN, TCMAX, ITC 


XDMIN, XDMAX, I XD 


WRITE(*,1100) 
READ (*,*) TL 


DH =(IZ-NRDOT)*(MASS-YVDOT)- 

& (MASS *XG-YRDOT) *( MASS *XG-NVDOT ) 
AA11=((IZ-NRDOT) *Y¥V-(MASS *XG-YRDOT) *NV)/DH 
AA12=((IZ-NRDOT) *(-MASS+YR)- 

& (MASS *XG-YRDOT) *(-MASS*XG+NR) )/DH 
AA21=((MASS-YVDOT) *NV-(MASS *XG-NVDOT) *YV)/DH 
AA22=((MASS-YVDOT) * (-MASS*XG+NR)- 

& (MASS *XG-NVDOT) *(-MASS+YR))/DH 
BB11=( (1Z-NRDOT) *YDRS—( MASS *XG-YRDOT) *NDRS ) /DH 
BB12=((IZ-NRDOT) *YDRB-(MASS*XG-YRDOT ) *NDRB ) /DH 
BB21=( (MASS-YVDOT) *NDRS- (MASS *XG-NVDOT) *YDRS ) /DH 
BB22=( (MASS-YVDOT) *NDRB-(MASS*XKG-NVDOT) *YDRB)/DH 


WRITE (*,1004) 


on 


10 


READS), 9 RATIO 


BB1=BB11+RATIO*BB12 
BB2=BB21+RATIO*BB22 


EPS =1.D-5 
ILMAX=1500 


DO 1 I=1,1TC 
WRITE (*,2001) I,1TC 
TC=TCMIN+ (1-1) *(TCMAX-TCMIN)/(ITC-1) 
TCD=TC*L/U 
POLEC=1.0/TCD 
AD1=3.0*POLEC 
AD2=3.0*POLEC**2 
AD3=POLEC** 3 
Al=BB1*U*U 
B1=BB2*U*U 
Cl=-AD1-(AA11+AA22)*U 
A2=(BB1*AA22-BB2*AA12) *U** 3 
B2=(BB2*AA11-BB1*AA21 ) *U**3 
K1=AD3/((BB2*AA11-BB1*AA21) *U**3) 
C2=AD2-(AA11*AA22-AA12*AA21) *U**2+BB2*U*U*K1 
K2=(C1*B2-C2*Bl)/(Al*B2-A2*B1 ) 
K3=(C2*A1-Cl*A2)/(A1*B2-A2*B1) 
DO 2 J=1,IXD 

XD=XDMIN+(J~1) *(XDMAX-XDMIN) /(1XD-1) 


CALL LINEAR(K1,K2,K3,U,XD,AA11,AA12,AA21,AA22,BB1,BB2,A,TL) 
CALL RG(4,4,A,WR,W1,0,222,1V1,FV1,IERR) 
CALL DSTABL(DEOS,WR,WI, FREQ) 


IF (J.GT.1) GO TO 10 
DEOSOO=DEOS 

XDOO =XD 

LL=0 

GO TO 2 

DEOSNN=DEOS 

XDNN =XD 
PR=DEOSNN*DEOSOO 

IF (PR.GT.0.D0) GO TO 3 
LL=LL+t+l 

IF (LL.GT.3) STOP 1000 
IL=0 

XDO=XDOO 

XDN=XDNN 
DEOSO=DEOSOO 
DEOSN=#=DEOSNN 

XDL=XDO 

ADR=XDN 

DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 


CALL LINEAR(K1,K2,K3,U,XD,AA11,AA12,AA21,AA22,BB1,B8B2,A,TL) 
CALL RG(4,4,A,WR,WI,0,222,1V1,FV1, 1ERR) 
CALL DSTABL(DEOS,WR,WI, FREQ) 


DEOSM=DEOS 
XDM=XD 
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PRL=DEOSL*DEOSM 
PRR=DEOSR*DEOSM 


IF (PRL.GT.0.D0) GO TO 5 


XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 
DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 
DIF=DABS (XDL-XDM) 

IF (DIF.GT.EPS) GO TO 
XD=XDM 

Go TO 4 

IF (PRR.GT.0.D0) STOP 
XDO=XDM 

XDN=XDR 

DEOSO=DEOSM 
DEOSN=DEOSR 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 
DIF=DABS(XDM-XDR) 

IF (DIF.GT.EPS) GO TO 
XD=XDM 

LLL=10+LL 

WRITE (LLL,*) XD/L,TC 
XDOO=XDNN 
DEOSOO=DEOSNN 


CONTINUE 
CONTINUE 


FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 


ENTER 
ENTER 
ENTER 
ENTER 
ENTER 


MIN, MAX, 
MIN, MAX, 
U') 


(! 
2 
(2 BOW/STERN 
(! 


FORMAT (215) 


END 


3100 
6 


AND INCREMENTS OF TC‘) 
AND INCREMENTS OF XD‘°) 


RUDDER RATIO’ ) 


TIME LAG TL’) 


SUBROUTINE DSTABL(DEOS,WR,WI,OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4),WI(4) 
DEOS=-1 .0D+20 
DO 1 I#1,4 
IF (WR(I).LT.DEOS) GO TO 1 
DEOS=WR(I) 
IJ=!I 
CONTINUE 
OMEGA=WI (IJ) 
OMEGA=DABS( OMEGA) 
RETURN 


END 


SUBROUTINE LINEAR(K1,K2,K3,U,XD,AA11,AA12,AA21,AA22,BB1,BB2,A,TL) 


IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
DOUBLE PRECISION K1,K2,K3 
DIMENSION A(4,4) 


A(l, 
,2)=0.0D0 
,3)=1.0D0 
,4)=0.0D0 


A(1 
A(1 
A(l 


1)=0.0D0 
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A(2,1)= BB1*U*U*K1* (XD-U*TL) /XD 
A(2,2)=U*(BB1*U*(K2*XD-K1*TL)+AA11*XD) /XD 
A(2,3)=AA12*U+BB1 *U*U*K3 

A(2,4)= BB1*U*U*K1/XD 

A(3,1)= BB2*U*U*K1*(XD-U*TL) /XD 
A(3,2)=U*(BB2*U*(K2*XD-K1*TL)+AA21*XD) /XD 
A(3,3)=AA22*U+BB2*U*U*K 3 

A(3,4)= BB2*U*U*K1/XD 

A(4,1)=U 

A(4,2)=1.0D0 

A(4,3)2#0.0D0 

A(4,4)=0.0D0 

RETURN 
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APPENDIX B 


PROGRAM THESIS2.FOR (Time Delay-2nd Order Approx) 
BIFURCATION ANALYSIS 


IMPLICIT DOUBLE PRECISION (A-H,O-Z) 

DOUBLE PRECISION K1,K2,K3,L,NR,NV,NDRS,NDRB,IZ,MASS, 
& NRDOT, NVDOT 

DIMENSION A(4,4),FV1(4),1V1(4),222(4,4),WR(4),WI(4) 


OPEN (11,FILE='BIF1.RES’ ,STATUS=’ NEW’ ) 
OPEN (12,FILE='BIF2.RES’ ,STATUS=‘ NEW’ ) 
OPEN (13, FILE='BIF3.RES’ ,STATUS=‘NEW’ ) 


Vehicle Parameters: 
WEIGHT=435.0 


1Z =45.0 

L «7.3 
RHO =1.94 

G =32 22 
XG =0.0104 


MASS =WEIGHT/G 


YRDOT =-0.00178*0.5*RHO*L**4 
YVDOT =-0.03430*0.5*RHO*L** 3 
YR =+0.01187*0.5*RHO*LA* 3 
YV =-0.03696*0.5*RHO*L**2 
YDRS =4+0.02345*0.5*RHO*L**2 
YDRB =+0.02345*0.5*RHO*L**2 
NRDOT =-0.00047*%0.5*RHOAL#45 
NVDOT =-0.00178*0.5*RHO*L**4 
NR =-0.01022*0.5*RHO*L**4 
NV =-0.00769*0.5*RHO*L** 3 
NDRS =-0.337*0.02345*0.5*RHO*L**3 
NDRB =+0.283*0.02345*0.5*RHO*L*&* 3 


WRITE (*,1001) 

READ (*,*) TCMIN, TCMAX, ITC 
WRITE (*,1002) 

READ (*,*) XDMIN, XDMAX, IXD 
XDMIN=XDMIN*L 

XDMAX=XDMAX *L 

WRITE (*,1003) 

READ (*,*) U 


WRITE(*,1100) 
READ (*,*) TL 


DH =(IZ-NRDOT)*(MASS-YVDOT)- 

& (MASS*XG-YRDOT) *(MASS*XG-NVDOT ) 
AA11=((1IZ-NRDOT) *YV-(MASS*XG-YRDOT) *NV) /DH 
AA12=((1Z-NRDOT) *({-MASS+YR)- 

& (MASS *XG-YRDOT) *(-MASS*XG+4+NR))/DH 
AA21=( (MASS-YVDOT) *NV-(MASS*XG-NVDOT) *YV) /DH 
AA22=((MASS-YVDOT) *(-MASS*XG+NR) - 

& (MASS *XG-NVDOT ) *(-MASS+YR))/DH 
BBl11l=((IZ-NRDOT) *YDRS-(MASS*XG-YRDOT) *NDRS) /DH 
BB12=((1Z-NRDOT) *YDRB-(MASS*XG-YRDOT) *NDRB) /DH 
BB21=((MASS-YVDOT ) *NDRS-(MASS*XG-NVDOT) *YDRS) /DH 
BB22=((MASS-YVDOT) *NDRB-(MASS*XG-NVDOT) * YDRB) /DH 


WRITE (*,1004) 


Wl 


10 


READ (*,*) RATIO 


BB1=BB11+RATIO*BB12 
BB2=BB21+RATIO*BB22 


EPS =1.D-5 
ILMAX=1500 


DO 1 I=1,ITC 
WRITE (*,2001) I,ITC 
TC=TCMIN+(I-1)*(TCMAX-TCMIN) /( ITC-1) 
TCD=TC*L/U 
POLEC#=1.0/TCD 
AD1=3.0*POLEC 
AD2=3.0*POLEC**2 
AD3=POLEC ** 3 
Al=BB1*U*U 
B1=BB2*U*U 
Cl=-AD1-(AA11+AA22)*U 
A2=(BB1*AA22-BB2*AA1 2) *U** 3 
B2=(BB2*AA11-BB1*AA21 ) *U**3 
K1L=AD3/((B8B2*AA11-BB1*AA21 ) *U**3) 
C2=AD2-(AA11*AA22-AA12*AA21 ) *U**24+BB2*U*UAK1 
K2=(C1*B2-C2*B1)/(A1*B2-A2*B1 ) 
K3=(C2*A1-C1*A2)/(A1*B2-A2*B1 ) 
DO 2 J=1,IXD 

XD=XDMIN+(J-1) *( XDMAX-XDMIN)/(IXD-1) 


CALL LINEAR(K1,K2,K3,U,XD,AA11,AA12,AA21,AA22,881,BB2,A,TL) 
CALL RG(4,4,A,WR,WI,0,222Z,1V1,FV1,1ERR) 
CALL DSTABL(DEOS,WR,WI, FREQ) 


IF (J3-GT.1) GOure 10 
DEOSOO=DEOS 

XDOO =XD 

LL=0 

GO TO 2 

DEOSNN=DEOS 

XDNN =XD 
PR=DEOSNN*DEOSOO 

IF (PR.GT.0.D0) GO TO 3 
LL=LL+1 

IF (LL.GT.3) STOP 1000 
IL=0 

XDO=XDOO 

XDN=XDNN 
DEOSO=DEOSOO 
DEOSN=DEOSNN 

XDL=XDO 

XDR=XDN 

DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 


CALL LINEAR(K1,K2,K3,U,XD,AA11,AA12,AA21,AA22,BB1,BB2,A,TL) 
CALL RG(4,4,A,WR,WI,0,Z222,1V1,FV1,IERR) 
CALL DSTABL(DEOS,WR,WI, FREQ) 


DEOSM=DEOS 
XDM=XD 


ee 


PRL=DEOSL*DEOSM 
PRR=DEOSR*DEOSM 


IF (PRL.GT.0.D0) GO TO 5 


XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 
DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 
DIF=DABS (XDL-XDM) 

IF (DIF.GT.EPS) GO TO 
XD=XDM 

GO TO 4 

IF (PRR.GT.90.D0) STOP 
XDO=XDM 

XDN=XDR 

DEOSO=DEOSM 
DEOSN=DEOSR 

IL=IL+1l 

IF (IL.GT.ILMAX) STOP 
DIF=DABS(XDM-XDR) 

IF (DIF.GT.EPS) GO TO 
XD=XDM 

LLL=10+LL 

WRITE (LLL,*) XD/L,TC 
XDOO=XDNN 
DEOSOO=DEOSNN 


CONTINUE 
CONTINUE 


FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 


ENTER 
ENTER 
ENTER 
ENTER 
ENTER 


MIN, 
MIN, 
U’) 

BOW/STERN 


MAX, 
MAX, 


~= ~_ » = © 


FORMAT (215) 


END 
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AND INCREMENTS OF TC’) 
AND INCREMENTS OF XD’ ) 


RUDDER RATIO’ ) 


TIME LAG TL’) 


SUBROUTINE DSTABL(DEOS,WR,WI,OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-2Z) 
DIMENSION WR(4),WI(4) 
DEOS=-1 .0D+20 
DO 1 I=1,4 
IF (WR(I).LT.DEOS) GO TO 1 
DEOS=WR(TI) 
IJ=I 
CONTINUE 
OMEGA=WI (IJ) 
OMEGA=DABS ( OMEGA) 
RETURN 


END 


SUBROUTINE LINEAR(K1,K2,K3,U,XD,AA11,AA12,AA21,AA22,BB1,BB2,A,TL) 


IMPLICIT DOUBLE PRECISION (A-H,0-2Z) 
DOUBLE PRECISION K1,K2,K3 
DIMENSION A(4,4) 


A(1, 
A(1, 
A(1, 
A(1, 


1)=0.0D0 
2)=0.0D0 
3)=1.0D0 
4)=0.0D0 


es 


A(2,1)=(BB1L*U*U*K1-BB1*U*U*U* KL *&TL/XD) / 


& (1-BBL*U*U*KLATL*TL/(2.D0*XD) ) 
A(2,2)=(AA11*U+BB1 *U*U*K2-BB1L*U*U*K1*TL/XD) / 

& (1-BB1*U*U*KL*TL*TL/(2.D0*XD) ) 
A(2,3)=(AA12*U+BB1*U*U* K3+BB1*U*U*U* KIATL*TL/(2.D0*XD) )/ 

& (1-BB1*U*U*KI*TLATL/(2.D0*XD) ) 
A(2,4)=(BB1*U*U* K1/XD) / 

& (1-BBL*U*UAKLATLATL/(2.D0*XD) ) 
A(3,1)=(BB2*U*U*K1-BB2*U*U*U*KI*TL/XD) + 

& (BB2*U*U*K1L*TL*TL/(2.D0*XD) )* 

& (BB1*U*U*K1-BB1L*U*U*U* KL *TL/XD) / 

& (1-BBL*U*U*KLATL*TL/(2.D0*XD) ) 
A(3,2)=(AA21*U+BB2*U*U*K2-BB2*U*U*K1*TL/XD) + 

& (BB2*U*U*AKL*TL*TL/(2.D0*XD) )* 

& (AA11*U+BB1*U*U*K2-BB1*U*U*K1*TL/XD) / 

& (1-BB1*U*U*K1*TL*TL/(2.D0*XD) ) 
A(3,3)=(AA22*U+BB2*U*U*K3+BB2*U*U*UAKL*TL*TL/(2.D0*XD) )+ 

& (BB2*U*U*K1*TL*TL/(2.D0*XD) ) * 

g (AA12*U+BB1 *U*U*K3+BB1 *U*U*U* KL *TL*TL/(2.D0*XD) )/ 


(1-BB1*U*U*K1*TL*TL/(2.D0*XD) ) 
A(3,4)=(BB2*U*U*K1/XD) +(BB2*U*U*K1L*TLATL/(2.D0*XD) ) * 
E (BBL*U*U*K1/XD)/ 


& (1-BB1*U*UAKI*TL*TL/(2.D0*XD) ) 
A(4,1)=2U 

A(4,2)=1.0D0 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

RETURN 

END 
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APPENDIX C 


PROGRAM THESIS3.FOR (Time Delay-3rd Order Approx) 
BIFURCATION ANALYSIS 


IMPLICIT DOUBLE PRECISION (A-H,0-2Z) 
DOUBLE PRECISION K1,K2,K3,L,NR,NV,NDRS,NDRB,1IZ2,MASS, 


& NRDOT,NVDOT 
DIMENSION A(5,5),B(5,5),FV1(5),1V1(5),222(5,5),ALFR(5), 
6 ALFI(5),BETA(5),WR(5),WI(5) 


OPEN (11,FILE='’BIF1.RES’ ,STATUS='NEW’ ) 
OPEN (12,FILE='BIF2.RES’ ,STATUS='NEW’ ) 
OPEN (13,FILE='BIF3.RES’ ,STATUS='NEW’ ) 


Vehicle Parameters: 
WEIGHT=435.0 


IZ =45.0 

L =7.3 
RHO =1.94 

G m32.2 
XG =0.0104 


MASS =WEIGHT/G 


YRDOT =-0.00178*0.5*RHO*L**4 
YVDOT =-0.03430*0.5*RHO*L** 3 
YR =+0.01187*0.5*RHO*L* *3 
YV =—0.03896*0.5*RHO*L**2 
YDRS #+0.02345*0.5*RHO*L**2 
YDRB #=+0.02345*0.5*RHO*L**2 
NRDOT =-0.00047*0.5*RHO*L**5 
NVDOT =-0.00178*0.5*RHO*L**4 
NR =—0 .01022*0.5*RHO*L**4 
NV =—0.00769*0.5*RHO*L**3 
NDRS ©=-0.337*0.02345*0.5*RHO*L** 3 
NDRB =+0.283*0.02345*0.5*RHO*L** 3 


WRITE (*,1001) 

READ (*,*) TCMIN, TCMAX, ITC 
WRITE (*,1002) 

READ (*,*) XDMIN, XDMAX, IXD 
XDMIN=XDMIN*L 

XDMAX=XDMAX*L 

WRITE (*,1003) 

READ (*,*) U 

WRITE(*,1100) 

READ (*,*) TL 


DH =(I2-NRDOT)*(MASS-YVDOT)- 

& (MASS *XG-YRDOT) * (MASS *XG-NVDOT) 
AA11=((1IZ2-NRDOT) *YV-(MASS*XG-YRDOT) *NV) /DH 
AA12=((IZ2—-NRDOT) *(-MASS+YR)- 

& (MASS *XG-YRDOT ) *(-MASS*XG+NR) ) /DH 
AA21=((MASS-YVDOT) *NV-(MASS*XG-NVDOT) *YV) /DH 
AA22=((MASS-YVDOT) *(-MASS*XG+NR)- 

& (MASS*XG-NVDOT) *(-MASS+YR) )/DH 
BB11=((1I2-NRDOT) *YDRS-(MASS*XG-YRDOT) *NDRS)/DH 
BB12=((IZ-NRDOT) * YDRB-(MASS*XG-YRDOT) *NDRB) /DH 
BB21=( (MASS-YVDOT ) *NDRS-(MASS*XG-NVDOT) *YDRS) /DH 
BB22=( (MASS-YVDOT) *NDRB-(MASS*XG-NVDOT) *YDRB) /DH 


WRITE (*,1004) 


dis 


READ (*,*) RATIO 


BB1=BB11+RATIO*BB12 
BB2=BB21+RATIO*BB22 


EPS =1.D-5 
ILMAX=1500 


DO 1 I=1,1TC 
WRITE (*,2001) I,ITC 
TC=TCMIN+(I-1)*(TCMAX-TCMIN) /(ITC-1) 
TCD=TC*L/U 
POLEC=1.0/TCD 
AD1=3.0*POLEC 
AD2=3.0*POLEC**2 
AD3=POLEC**3 
Al=BB1*U*U 
B1l=BB2*U*U 
Cl=-AD1-(AA11+AA22)*U 
A2=(8B1*AA22-BB2*AA12) *U**3 
B2=(BB2*AA11-BB1*AA21) *U**3 
K1=AD3/((BB2*AA11-BB1*AA21) *U**3) 
C2=AD2-(AA11*AA22-AA12*AA21 ) *U**24+BB2*U*U*K1 
K2=(C1*B2-C2*B1)/(Al1*B2-A2*B1) 
K3=(C2*A1-C1*A2)/(A1*B2-A2*B1 ) 
DO 2 J=1,I1XD 

XD=XDMIN+(J-1) *(XDMAX=XDMIN) /( IXD-1) 


CALL LINEAR(K1,K2,K3,U,XD,AA11,AA12,AA21,AA22,BB1,BB2,A,B,TL) 
CALL RGG(5,5,A,B,ALFR,ALFI,BETA,0,222, TERR) 
DO 11 IJE=1,5 
WR(IJE)=ALFR(IJE)/BETA(IJE) 
WI(IJE)=ALFI(IJE)/BETA(IJE) 
CONTINUE 
CALL DSTABL(DEOS,WR,WI, FREQ) 


IF (J.GT.1) GO TO 10 
DEOSOO=DEOS 

XDOO =XD 

LL=0 

GO TO 2 

DEOSNN=DEOS 

XDNN «&=XD 

PR=DEOSNN* DEOSOO 

IF (PR.GT.0.D0) GO TO 3 
LL=LL+1l 

IF (LL.GT.3) STOP 1000 
IL=0 

XDO=XDOO 

XDN=XDNN 
DEOSO=DEOSOO 
DEOSN=DEOSNN 

XDL=XDO 

XDR=XDN 

DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 


CALL LINEAR(K1,K2,K3,U,XD,AA11,AA12,AA21,AA22,BB1,BB2,A,B,TL) 
CALL RGG(5,5,A,B,ALFR,ALFI,BETA,0,2Z22Z, IERR) 
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iZ 


DO 12 IJE1,5 


WR( IJE)=ALFR(IJE)/BETA(IJE) 
WI(IJE)=ALFI(IJE)/BETA(IJE) 


CONTINUE 


CALL DSTABL(DEOS,WR,WI, FREQ) 


DEOSM=DEOS 
XDM=XD 
PRL=DEOSL*DEOSM 
PRR=DEOSR*DEOSM 


IF (PRL.GT.0.D0) GO TO 5 


XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 
DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 
DIF=DABS(XDL-XDM) 

IF (DIF.GT.EPS) GO TO 
XD=XDM 

GO TO 4 

IF (PRR.GT.0.D0) STOP 
XDO=XDM 

XDN=XDR 

DEOSO=DEOSM 
DEOSN=DEOSR 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 
DIF=DABS( XDM-XDR) 

IF (DIF.GT.EPS) GO TO 
XD=XDM 

LLL=10+LL 

WRITE (LLL,*) XD/L,TC 
XDOO=XDNN 
DEOSOO=#DEOSNN 


CONTINUE 
CONTINUE 


3100 


FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT (’ ENTER TIME 
FORMAT (215) 

END 


ENTER MIN, 
ENTER MIN, 
ENTER U’ ) 


=-_ * -_ 


oe! gE gr or 


MAX, AND INCREMENTS OF TC’) 
MAX, AND INCREMENTS OF XD’) 


ENTER BOW/STERN RUDDER RATIO’ ) 


LAG TL’ ) 


SUBROUTINE DSTABL(DEOS,WR,WI,OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0O-Z) 


DIMENSION WR(5),WI(5) 
DEOS=-1.0D+20 
DO 1 I=1,5 
IF (WR(I).LT.DEOS) 
DEOS=WR(I) 
IJ=I 
CONTINUE 
OMEGA=WI (IJ) 
OMEGA=DABS( OMEGA) 
RETURN 
END 


GO TO l 


Wy 


SUBROUTINE LINEAR(K1,K2,K3,U,XD,AA11,AA12,AA21,AA22,BB1,BB2, 
& A,B,TL) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION K1,K2,K3 

DIMENSION A(5,5),B(5,5) 
A(1,1)=+0.0D0 
0.0D0 
1.0D0 
0.0D0 
0.0D0 
0.0D0 
0.0D0 
0.0D0 
0.0D0 
1.0D0 
(BB1*U*U*K1-BB1*U*U*U*K1*TL/XD) 
(AA1L11*U+BB1*U*U*K2-BB1L*U*U*K1*TL/XD) 
( 

( 

( 
U 

1 

0 

0 

0 

( 

( 

( 

( 

( 


» 
oe, 
— 
= 


AA12*U+BB1L*U*U*K3+BB1*U*U*AU*A KIA TLATL/(2.D0*XD) ) 
BB1*U*U*K1/XD) 
-1.D0+BB1*U*U*KI*TL*TL/(2.D0*XD) ) 


.0D0 

.O0D0 

.0D0 

.0D0 

BB2*U*U*K1-BB2*U*U*U*KL*TL/XD) 
AA21*U+BB2*U*U*K2-BB2*U*U*K1I*TL/XD) 
AA22*U+BB2*U*U*K3+BB2*U*UAUAKIATL*TL/(2.D0*XD) ) 
BB2*U*U*K1/XD) 

BB2*UAU*KI*TLATL/(2.D0*XD) ) 


MMrnNnNnNnh bb bh & Ww WW Ww WW NM NN DN FF Fe 
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.0D0 
.O0D0 
.0D0 
.0D0 
.0D0 
.0D0 
.ODO 
.0D0 
.0D0 
.OD0 
.0D0 
.0D0 
BBL*U*U*U* KIA TLATLATL/(6.D0*XD) ) 
.O0D0 
BBL*U*UAKL*TL*TLATL/(6. DO*XD) ) 
.0D0 
.0D0 
.0D0 
.0D0 
.0D0 
.0D0 
.0D0 
1.0D0+BB2*U*U*U*KI*TLATL*TL/(6.D0*XD) ) 
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B(5,5)=(BB2*U4U*K1*TL*TL*TL/(6.D0*XD) ) 
RETURN 
END 
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APPENDIX D 


PROGRAM THESISN.FOR (Time Delay - Newton’s Iteration Method) 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION K1,K2,K3,L,NR,NV,NDRS,NDRB,1IZ,MASS, 
& NRDOT,NVDOT 

DIMENSION A(4,4),FV1(4),1V1(4),222(4,4),WR(4),WI(4) 


OPEN (11,FILE=’BIF1.RES’ ,STATUS=’ NEW’ ) 


Vehicle Parameters: 
WEIGHT=435.0 


1Z =45.0 

L =] .3 
RHO =1.94 

G =32.2 
XG =0.0104 


MASS =WEIGHT/G 


YRDOT =-0.00178*0.5*RHO*L**4 
YVDOT =-0.03430*0.5*RHO*L** 3 
YR =+0.01187*0.5*RHO*AL** 3 
X¥V =-0.03896*0.5*RHO*L*42 
YDRS =+0.02345*0.5*RHO*L**2 
YDRB =+0.02345*0.5*RHO*L**2 
NRDOT =-0.00047*0.5*RHO*L**5 
NVDOT =-0.00178*0.5*RHO*L**4 
NR =-0.01022*0.5*RHO*L**4 
NV =-0.00769*0.5*RHO*L**3 
NDRS =-0.337*0.02345*0.5*RHO*L«*3 
NDRB =+0.283*0.02345*0.5*RHO*L** 3 


WRITE (*,1001) 


READ (*,*) TCMIN, TCMAX, ITC 
WRITE (*,1003) 
READ (*,*) U 


DH =(1Z-NRDOT)*(MASS-YVDOT)- 

& (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 
AA11=((1Z-NRDOT) *YV-(MASS*XG-~YRDOT) *NV)/DH 
AA12=((1Z-NRDOT) *(-MASS+YR)-~ 

& (MASS*XG-YRDOT) * (-MASS*XG+NR) ) /DH 
AA21=((MASS-YVDOT) *NV-(MASS*XG-NVDOT) *YV) /DH 
AA22=( (MASS-YVDOT) *(-MASS*XG+NR) - 

& (MASS*XG-NVDOT) *(-MASS+YR) )/DH 
BB11=( ( IZ-NRDOT) *YDRS~(MASS*XG-YRDOT) *NDRS) /DH 
BB12=((1Z-NRDOT) *YDRB-(MASS*XG-YRDOT ) *NDRB) /DH 
BB21=((MASS-YVDOT) *NDRS-(MASS*XG-NVDOT) *YDRS) /DH 
BB22=( (MASS-YVDOT) *NDRB-(MASS*XG-NVDOT) *YDRB) /DH 


WRITE (*,1004) 
READ (*,*) RATIO 


BB1=BB11+RATIO*BB12 
BB2=BB21+RATIO*BB22 


WRITE (*,1005) 
READ (*,*) TL 


EPS=1.D-10 
ITMAX=1000 


(2 


AAD 


aan 


Rm BM mM 


DOV = le 


WRITE (*,2001) I,1TC 
TC=TCMIN+(I-1)*( TCMAX-TCMIN)/(ITC-1) 
TCD=TC*L/U 

POLEC#1 .0/TCD 

AD1=3.0*POLEC 

AD2=3.0* POLEC**2 

AD3=POLEC®% * 3 

A1=BB1L*U*U 

Bl=BB2*U*U 

Cl=-AD1-(AA11+AA22) *U 
A2=(BB1*AA22-BB2*AAI 2) *U** 3 
B2=(BB2*AA11-BB1*AA21 ) *UX* 3 
K1=AD3/((BB2*AA1L1-BB1*AA21 )*U**3) 
C2=AD2-(AA1L1*AA22-AA12*AA21 ) *U**2+BB2*U*U* KI 
K2=(C1*B2-C2*B1)/(A1L*B2-A2*B1 ) 
K3=(C2*A1-C1L*A2)/(A1L*B2-A2*B1 ) 


Al= (AA11*BB2-AA21*BB1 )*U*U*U* KI 

A2= (AA1L1L*AA22-AA12*AA21 )*U*U+(AA11*BB2-AA21*BB1) *U*U*U*K 3+ 
(AA22*BB1-AA12*BB2) *U*U*U* K2-BB2*U*U* KI 

A3=-(AAL1+4+AA22)*U-(BB1*K2+BB2*K3) *U*U 

BO= (AA11*BB2-AA21*BB1 ) *U*U*U*U4 KI 

Bl=-—(AA12*BB2-AA22*BB1) *U*U*U* KIL-BB2*U*U*eUe KL 

B2=-BB1*U*U* KI 


COMPUTE INITIAL APPROXIMATION FOR OMEGA (TL#=0) 


AQ=B2*A3-Bl 
BQ=B1*A2-B2*A1-B0*A3 

CQ=BO*AlL 

DET=BQ* *2-—4.D0*AQ*CQ 

IF (DET LT.0.DO0O)° GO Ton: 
ZT1=(-BQ+DSQRT( DET) )/(2.0D0*AQ) 
ZT2=(-BQ-DSQRT(DET))/(2.0D0*AQ) 
IF (ZTL.GE.0.D0) OM=DSQRT(ZT1) 
IF (ZT2.GE.0.D0) OM=DSQRT(ZT2) 
ALPHA1L=AQ 

ALPHA2=BQ 

ALPHA3=CQ 

BETA1=-B2 

BETA2=B0+B2*A2-B1*A3 
BETA3=B1*A1-B0*A2 


COMPUTE EXACT OMEGA USING NEWTON'S ITERATION 


OMOLD=0M 

F=(BETA1 *OMOLD* *5+BETA2 *OMOLD®* * 3+BETA3* OMOLD ) *DSIN( OMOLD* TL) 
+(ALPHA1*OMOLD®* * 4+ALPHA2*OMOLD* * 2+ALPHA3 ) *DCOS( OMOLD* TL) 

FPRIME=(5.D0*BETA1 *OMOLD** 4+3.D0*BETA2*OMOLD* * 2+BETA3 ) 
*DSIN(OMOLD*TL)+(BETA1 *OMOLD* *5+BETA2 *OMOLD* *3+BETA3*OMOLD) 
*TL*DCOS(OMOLD*TL)+(4.D0* ALPHA] *OMOLD**3+4+2.D0*ALPHA2*OMOLD ) 
*DCOS (OMOLD* TL) —-( ALPHA1 *OMOLD* * 4+ALPHA2 *OMOLD* *2+ALPHA 3 ) 
*TL*DSIN( OMOLD*TL) 

IF (FPRIME.EQ.0.00) STOP 1112 

IF (F.EQ.0.D0) GO TO 2 

OMNEW=OMOLD-F/FPRIME 

OMDIF=DABS ( OMNEW-OMOLD ) 

IF (OMDIF.LE.EPS) GO TO 2 
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OMOLD=OMNEW 

IT=IT+1 

IF (IT.GT.ITMAX) STOP 1111 
GO TO 3 


OM=OMNEW 

XDNUM=DSQRT( (BO0-B2*OM**2) **2+(B1*OM)**2) 
XDDEN=OM*DSQRT( (A1-A3*OM**2 )**2+4(A2*OM-OM** 3) **2) 
XD=XDNUM/XDDEN 

XD=XD/L 


WRITE(11,*) XD,TC,OM 
CONTINUE 


FORMAT 
FORMAT 


(‘ ENTER MIN, MAX, AND INCREMENTS OF TC’) 

( 
FORMAT ( 

( 

( 


t 

‘ ENTER U‘) 

‘ ENTER BOW/STERN RUDDER RATIO’ ) 
FORMAT (‘ ENTER TIME LAG‘ ) 
FORMAT (215) 


END 
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